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Introduction 


I„„a St.ee University Engineering Research Institute Project 9U-S is 

a continuation of RASA Grant BGRIG-002-029 . The major area, of study and 

goals of this research project are: 

„ Determine the flow field generated by a finitely thick lilting 

three-dimensional wing with subsonic tips moving at superson.c 
speeds. 

2) Determine the flow field produced by a lilting elliptic cone 

using finite difference techniques. 

3) study the problem of cross-flow instability associated with lift- 
ing delta wing configurations such as space shuttle. 

A sis month extension of the current grant through December 1972 has 
been obtained in order to complete the first objective of the research project. 
Result, obtained during the past year are presented through theses or reports 

where applicable. 


1. Rectangular Wing Problem 


The flow field associated with the forebody region of the symmetric 
rectangular wing has been determined for angles-of-attack of 0 and 4 
degrees at a Mach number of 2. The original solution technique has been 
modified to include the entropy boundary conditions developed by Mike 
Abbett, thus replacing the reflection technique in the evaluation of 
the ultimate body values of the flow variables. Although the application 
of the boundary conditions as described by Abbett failed to yield a 
viable solution in the cone region (the wedge region solution was satisfactory) 
a suitable modification was made that did work. It would appear that in 
regions of small radius of curvature (in the cone-tip region) the predicted 
values of the velocity components on the body are not sufficiently accurate 
to use as inputs for the Abbett technique which would continuously expand 
the flow on the body. The modification consisted of using corrected body 
data obtained via reflection as inputs to the Abbett technique. 

The flow field associated with the afterbody region of the rectangular 
wind is currently being developed. The forebody flow field data is transferred 
to the afterbody grid system at the wind midchord point, thus generating 
initial data for the afterbody. The integration technique and boundary 
conditions used are the same as for the forebody. To date, the integration 
has proceeded to a point approximately 90% of the distance back on the 
afterbody. In this region of the flow field on the body the cross flow 
eigenvalues grow rapidly as the trailing edge is encountered forcing the 
integration step size to be continuously decreased. Since the normal 
eigenvalues do not change appreciably along the afterbody the effective 
Courant number at which the solution is being obtained in the normal 
direction continuously decreases yielding a suboptimal solution in the 
normal direction in the region near the trailing edge of the win&. The 
situation has been improved somewhat by decreasing the number of grid 
points about the cone tip. It would appear that the afterbody solutions 
will be completed in the very near future. 


The afterbody solutions will be used as initial data for the wake 
region integration. No problems are anticipated in the wake integration 
which should be completed shortly after the afterbody integration is 
furnished. 


2) Elliptic Cone Flow Fields 

Appendix A includes the thesis "Finit difference solutions to the 
equations of fluid flow around an elliptic cone" by Phil Pratt. This seiv- 
ed as Phil's M.S. thesis and deals with the second part of the current effort. 

The results presented compare favorably with both experimental and ana- 
lytical data available. The only question arises when one notes the persistent 
oscillation in body pressure which occurs approximately 115° around the ellip- 
tic cones both at zero and non-zero angles of attack. This question should 
have been answered by a more complete investigation. However, the funds avail- 
able for computer time and the time requirements precluded further investigation. 

Since Pratt's work was completed, a new program has been written for solving 
the elliptic cone flow field. This has been done at the direction of the NASA 
Technical Monitor Dr. Kutler and will be used both in further studies of the 
elliptic cone and the cross-flow instability problem. This program is current 
in the de-bug stages. A complete description of the new program will be presenter 
in a later report when computed results are available. 

3) Cross-Flcw Instability Problem 

The third part of the current effort is to study the problem of cross flow 
instability associated with very small radius of curvature areas such as the lead- 
ing edge of the shuttle wing. This phase of the study will begin after the com- 
pletion of the elliptic cone flew field study. The programs resulting from the 
elliptic cone study are necessary to the cross-flow investigation. 

A) Other Work Partially or Wholly Supported Under the Current nr ant 

A) Numerical Advection Experiments 
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Appendix B is a paper which has been submitted to the Journal of 
Applied Meteorology for publication. This work was partially funded 
through the current grant and is included here. Behrooz Fattahi was the 
Graduate Assistant working on that particular paper. 

B) Undefined Problems 

Jim Daywitt has been supported by th<* current grant during the past 
year and will be starting on his thesis in September. His thesis topic 
is still undefined. It will be in an appropriate area dealing with finite 
difference methods in aerodynamics. 
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A,B,C 


E,F,G, H 


AAA 

i, j.k 

A A A 

i„ , i , i 
R N C 


S 

u, v, v 
u|,^ ,w$ 
u 

x,y,z 


Coefficient matrices of the gas dynamic equations 

Semimajor axis of the body 

Semiminor axis of the body 

Characteristic length of the body 

Matrices of conservative variables ( Appendix A ) 

Identity matrix 

Unit vectors on the x,y,z directions 

Unit vectors in the radial, normal and cross flow 
directions 

Total enthalpy 

Mach number 

Dimensionless pressure 
Dimensionless velocity vector 

Dimensionless velocity vector in the g, direction 
Radial vector ( Piqure 1 ) 

Radius of circular cross section in computational 
system 

Entropy 

Dimensionless velocity components of q 

Dimensionless velocity components of cfj 

Vector of independent flow variables 

Coordinate directions in the physical Cartesian 
system 

Coordinate directions in the computational 
Cartesian system 
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0! 

Angle of attack 
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Ratio of specific heats 

Ap 

Increment, in the p-direction 

A6 

Increment in the 6- direction 

Af 

Increment in the F-direction 

AS 

Mesh ratio 

A6 


A? 

A3 

Mesh ratio 

A9 

Angle through which flow is 
( Appendix C ) 

6,3 

Coordinate directions in the 
system 

P 

Dimensionless density 

o, X 

Eigenvalues 

7 

Vector operator 


expanded or compresse 
computational polar 


d 


Subscr ipts: 


c 

Pertaining 

to the circular cone 

c 

Cress flow 

coordinate direction ( Figure 1 ) 

e 

Pertaining 

to the elliptic cone 

1 

Mesh point. 

in the 6-direction 

k 

Mesh point 

in the 3-direction 

max 

Maximum 


N 

Coordinate 
( Figure 1 

direction normal to the body surface 
) 


tan 


Radial coordinate direction ( Figure 1 ) 
Tangential component 


,z 


Partial derivatives with respect to x 
x ’ » Y ',z' Partial derivatives with respect to x 

T Flow variables after the predictor 

2 Flow vari able s after t he cor r ac t. o r 

oo Freestream condition 

p Partial derivative with respect to %, 

Superscripts: 

~ Predicted value cf flow variables or 

matrix forms 

n Spacial step location in integration 

♦ SuperlayeL, sublayer 

— Dimensional quantity 


, y or z 
• , y * or 7 . * 


6 or 3 

intermediate 


SUMMARY 


The solutions to the equations of notion for inviscid 
fluid flow around a pointed elliptic cone at incidence are 
presented. The numerical method used, MacCormack's second 
order preferential predictor-corrector finite difference ap- 
proximation, is applied to the fluid flow equations derived 
in conservation-law form. The entropy boundary condition, 
hitherto unused for elliptic cone problems, is investigated 
and compared to reflection boundary condition solutions. Th 
stagnation streamline movement of the inclined elliptic cone 
is noted and surface pressure coefficients are plotted. 

Also presented are solutions for an elliptic cone and a 
circular cone at zero incidence and a circular cone at a 
small angle of attack. Comparisons are made between these 
present solutions and previously published theory. 




introduction 




The reduction of complex bodies to simpler configura- 
tions in investigating aerodynamic chars -ter istics has long 
been accepted practice. Leading edges of thin airfoils have 
been assumed to be wedges of finite ( or infinite ) span, and 
fuselages have been assumed to be circular cones attached to 
cylinders. These assumptions enabled the sc ianti st /eng i neer 
to apply analytical methods to acquire solutions to the par- 
tial differential equations of fluid motion. However, as the 
aerodynamic "state of the art” grew, the configurations be- 
came more complex, and the '•simplifying” bodies became more 
complicated. Therefore, it was necessary to seek improved 
and more efficient methods of solution. 

The applicability of numerical methods to the solution 
of differential equations had been long known and understood, 
but it wasn't until the advent of the high-speed digital com- 
puter that efficient use could he made of these methods. Al- 
though these numerical approaches are quite amenable to many 
problem types, they should not be regarded as panaceas. Many 
times, more problems are created than solved with the intro- 
duction of the numerical techniques. Difficulties with nu- 
merical stability, accuracy, computer time and storage re- 
quirements, and applicability of the numerical metnods to the 
physical situation often cause serious problems. 
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This paper is primarily concerned with the numerical so- 
lution to nonlinear conical flow and, in particular, the flow 
around elliptic cones at incidence. The circular cone has 
been the accepted configuration for modeling the supersonic 
flow about bodies of revolution for decades. However, for 
non-axisymae t ric bodies immersed in a supersonic stream, the 
circular cone is of little use. Therefore, as Van Dyke ( 25 ) 
points out, the elliptic cone may become the model for these 
geometries. Furthermore, elliptic cones provide certain 
aerodynamic characteristics that surpass tnosa of thoir cir- 
cular counterpart { e.g., from Chapkis (7), greater L/D ratio 
for the same cross sectional area per unit height ). Evi- 
dence of the impact of the elliptic geometry can easily be 
seen in present-day lifting body and space shuttle designs. 

The problem or the supersonic flow field surrounding 
these bodies has been investigated by Perri (9,9,10) , Van 
Dyke (25), Helnick (19), Kahane and Solarski (13), Jones 
(12), Chapkis (7), and Martellucci (18), and others, with 
varying degrees of success. The papers by Fecri (10), 
Martellucci, Chapkis and Jones were of greatest assistance in 
this study, since they not only presented theoretical compu- 
tations, but also numerical and experimental comparisons. 

The publications by Ferri (d,9). Van Dyke, and Xahane and 
Solarski were basically of theoretical content. Mel nick's 
paper was of special interest, because it succinct ly ex- 
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plained the streamline pattern and vortical singularity prob- 
lems near the body. 

South and Klunker (24) indicate two general approaches 
to the solution of nonlinear conical flows. The first in- 
cludes the •'distance asymptotic" methods, where initial con- 
ditions are chosen near the conical apex, and numerical inte- 
gration proceeds downstream until conical similarity condi- 
tions ( defined later ) are satisfied. The second approach 
consists of "conical" methods that invoke the conical self- 
similarity, then reduce the problem to two independent vari- 
ables. There are different variations to the approaches 
above: Rakich (22), Chapkis, Martellucci and Ferri utilize 

different forms of the method of characteristics; Kahane and 
Solarski propose Ward's elliptic slender body theory; and 
South and Klunker implement the method of lines ( method of 
integral relations ). For a brief history of these develop- 
ments, the reader is referred to South and Klunker (24). 

The method chosen for this study was the "distance as- 
ymptotic" method, utilizing a finite difference approach. 
However, this approach was modified slightly to enhance sta- 
bility consider at ion s. Instead of integrating lownst.ream 
until convergence was attained, the integration was initiated 
at some point f , and a solution at sose point § ♦ was 
obtained. This solution was then used as the initial condi- 
tion for another integration step at point £ . A mote de- 
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tailed explanation of tbe reasoning and procedure is present- 
ed later. 

The elliptic body surface was defined mathematically, 
and boundary conditions were applied at the body surface and 
the unknown shock position. Freestream values for pressure, 
density, and velocities were used as initial conditions so 
that the numerical r epresentat ion of the shock build-up was 
as if the cone had been impulsively started. However, the 
transient numerical flow field behind the shock may not rep- 
resent the actual flow field until convergence has been 
attained. 

The equations of fluid motion were simplified by as- 
suming: 1) negligible viscosity, therefore eliminating 

boundary layer effects; 2) negligible heat conductivity ; 3) a 
perfect gas; and 4) conical flow, resulting in isentropic 
flow along streamlines. The physical properties of conical 
flow are such that pressure, density, velocities, energy and 
entropy remain constant along rays emanating from the conical 
apex ( i.e., the aforementioned conical similarity ) . 
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Co ordinate System 


The coordinate systea chosen, shown in Figure 1, was t<. 
result of a double coordinate transformation. The first 
transformation, from the physical Cartesian system, (x,y,z), 
to the computational Cartesian system, (x',y,z'|. converted 
the elliptic cross section to a circular cross section ( see 
Figure A- 1 ). This constituted a stretching oc the seinimajor 
and semiminor axes, a and b, of the elliptic cone by a factor 
of b and a, respectively; provided a,b > 0 . 

As a result of this transformation, angles around the 
body in the y,z-plane were not conserved in the y',z'-plane. 
Evenly spaced angles around the computat iona 1 Cartesian body 
produced rays that "clustered" about the 0 - } 0 o meridian in 
the physical flane ( see Figure A-1 ). The relation 
governing the angles is given by: 


where 0 c is the computational plane body angle. This behav- 
ior is beneficial for elliptic cones of Urge eccentricities. 
For these cases, flow near the 3 * 90« meridian is of great 
importance due to the presence of large gradients in the flow 
variables. Using this t ransforaation aLlows morv mesh points 
to be concentrated in this area, therefore monitoring the 
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flow nor* closely. 

The second transformation carried the computat ion al Cc" 
tesian plane into the computational polar plane, ( ?, 6 ) , 

seen in Figure 1. This system allowed for the easy applies 
tion of the elliptic boundary conditions to the circular v.one 
( described in more detail later ), while retaining the coni- 
cal similarity of the flow field. For the numerical integra- 
tion method chosen, requiring a coordinate axis to be coinci- 
dent with the conical axis ( i.e., x = x * - I ) yi^ldel 

simple functions of the gas dynamic variables in conservative 
form ( see Equation 1-11 >. These functions could then be 
easily solved to yield explicit functions of the gas dynamic 

variables. 

When the aforementioned assumptions ( 1) negligible 
Yiscosity; 2) negligible heat conductivity: 3) perfect gas; 

and «> conical flo. ) were applied to the equations of notion 
of fluid flo., they yielded the folloving equations for con- 

sideration: 


X-HomentUk : 

(P ♦ pu 2 )* ♦ ( puv) y ♦ ( Puw) z 

= 0 

(2) 

Y-nomentum: 

( puv) x ♦ (P ♦ P» l )» 4 ( pvw) z 

= 0 

(3) 

Z-Hoaentum: 

( pun) „ ♦ (pvw) Y ♦ (P ♦ ow 1 ) x 

= 0 

[<*) 


.7 # r op - c*X7o) 

= 0 

( c ) 


Energy: 


8 


Continuity: V*(Pq) - 0 . (M 

However, t horo equations ar* not yet in i f on that will 
allow numerical computation. They must, first h-' 1 i a nsf or me 1 
to ♦hn computarion.il polar cool"! incite system to b 1 : on si ; t on ♦ 

with t. ho geometry or t h * problon. Then, sin cc the tiinsfor- 
ib a t i on q i v e s r i ie to n o n h o m o q e n ? o u s terns, t. h => aquations must 
b»i expanse! to ronseivat ion-law form to accomodate the numer- 
ical method. This is done in Appendix For njm^rical sig- 

nificance, and as a matter of qeneral practice, it is conve- 
nient to nond imens i o na 1 i 2° t he flow parameters. The noni i- 
mensional izat ioii scheme u.se \ is presentM in Appendix c. 

N umer ical Method 

Th<-» numerical met ho! chosen was KacCorm ac k • ^ (16) second 

order preferential pre i ictor-corrector, utilizing a forward 
differencing predictor step and a backward differencing cor- 
rector Th i > method was .shewn by Kutler ( Q 9) to b- 5 su- 

perior to fiLst order an! many second order methods in terms 
ot computer time an! storage, programming * \ se and iccutacy. 
The method is i finite lilfereoc * techniju^ th it is terms! 

••;>: et * rent i il " dec t <i so t h- dif i •renc in g lirection in th* pi* 

.! i c t oi and rot i ■mM.ut > t * p . /i b ^ changed it will, provid'd 
t h«» .»* •:luu: jes ai t r.nieisfent with the second or 1 u iunuary, 
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When applied to a system of equations of form: 
E x ♦ F y + G z = 0 , 

HacCormack • s predictor-corrector is written as: 


V 

Predictor: 

e" + * „ 

E n - ~-(F n - F n ) 

J,k Ay j+l,k j,k 

Ax , 

* ■ I v-» 

A* j,k+l 

& 

f 


J,k 

V 

Corrector: 

E 0+1 = 

hi** ♦ F?*- ^(^ +I 

J,k i,k Ay j,k 

Ax ~n+l ~ n+ ’ 

- > 

‘v. ; 



i - 1 , k ’ 

• 

-v- 

1. 

s 

1 

However, as 

car. be 

seen from Appendix A, 

t h system 


tions has the form: 


E x ♦ Fy ♦ G z ♦ H = o 


To accomodate the nonhomoqenoous term, it was necessary to 
add the quantities: 


Predictor: - H n Ax 

^ j , k 

Corrector: - v7J n Ax 

j,k 


to Equations 8. 

For two-dimensional flow around circular cones, amplifi- 
cation theory predicts the numerical method is stable, pto 




viled ( uninq present notation ): 
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tantamount to tne satisfaction of Equation 12 . 
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This is ' a • to roe fact t. hat 
' 1 ■ '•>’ v i. ■ ••■as*.*d on linear F >urier analv 

• > t ;■ i -,i n oft i near . 

‘ * t ha t • ho < ; \ no initial • m 1 i t i or. ; a t 

' " ’ ii i» in i * > s ;oirt r * aC , and ♦ h-*n 

* : ' * • *''*' ; < ’1 n * i >r ■ i i ; 1 ' i i 1 r j n 1 1 1 i ) r.. . f o r 
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numerical 

accuracy. 

In effect, the numerical method 

was con- 

tinuously 

operating 

at ( or near ) the optimum mesh 

ratio. 

As can be 

seen from 

Appendix D: 



| a | , 

max 

1 ‘.ax' a — > 

(14) 






so that after a number of downstream integrations, the numer- 
ical method would no lonqer be functioning optimally. 

Initial Conditions 

Initial conditions for the flow variables can be chosen 
in at least two ways. First, the shock shape and position, 
and values of flew variables between the shock and the body 
can be approximated. This method is quite efficient with 
regard to computer effort, if the approximations are reason- 
able. Noretti (21) indicates that the choice of initial data 
does not affect, the end results. Secondly, the individual 
flow variables can be initialized to frsestream values. Al- 
though quite inefficient, this method was chosen since it 
resembles impulsively starting the body. Even though the 
transient, numerical values for the flow variables between 
the shock and body surface may not coincide with the true, 
physical solution, the shock build-up can still be witnessed. 
Appendix £ shows ho < the freestream quantities were derived. 
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Boundary Conditions 

From analytical considerations, the exact boundary con- 
dition applied to the surface of the body is: 

g • V F(?,6,P) = 0 , (15) 

where F (f.,6,0) 0 i; the ^qaation of th* 1 holy. Also, 

across the shock surface, th<-» Pankine-Hu.jor.iot equations must 
be satisfied. However, numerical method; cannot use these 
equations directly. This is especially true of finite dif- 
ference methods, du a to their 'differencing natur?. There- 
fore, equivalent numerical boundary conditions ware derived 
and applied at the body and shock. 

The reflection boundary condition jsed in this study 
( see Appendix D ), assumed even functions of pressure, den- 
sity, and tangential velocity and an odd function of normal 
velocity across tne bo dy surface. This gave t h a necessary 
sublayer data point.*;, as a function of the known superlayer 
points, requited by the KacCormack predictor ot-»p. Since el- 
liptic cones have a plan® of flow symmetry about the & - 0°, 
180° meridian ( sae Figure 1 ), numerical integration was 
cai tied out over one half the conical surface. "his plane of 
symmetry forms another boundary across which pressure, der.si- 
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ty, u-velocity, and v-velocity were assumed to be even func- 
tions and w-velocity was assumed to be an odd function. 

The mesh in the direction outward from the body, in the 
direction of increasing 6, extended far enough to overlap the 
estimated shock position. In this manner, the shock could, 
at best, be captured between two mesh points. Pressure, den- 
sity and velocity outside the shock were assumed to have 
freestream values. 

Due to first order numerical accuracy in the MacCormack 
predictor, the velocity vector, q, , was not, in general, tan- 
gent to the body surface as dictated by Equation 15. The en- 
tropy boundary condition ( see Appendix C ) used by Abbett 
(1) accounted for this misalignment and "brought” the veloci- 
ty vector back to tangency by utilizing a local compression 
or expansion through an angle A0. With A0 known, the 

pressure, P , following the expansion or compression was cal- 
2 

culated; and with pressure known, density, , at the body 
was obtained ty utilizing the perfect gas equation of state. 
From energy considerations, the velocity magnitude, | q | , was 
found. Finally, from vector relations inherent in the flow 
geometry, the direction of the resultant velocity, q , was 
obtained. A complete derivation is presented in Appendix c. 

For this problem, it was necessary to apply the entropy 
boundary condition after the predictor step. Mo analytical 
proof for this conclusion is readily available. However, it 


IS 


pel ? that the i mm e 1 ia t e application of th? entropy 




boundary condition ( i.e., after the praiictor ) allow; the 
con-'cti i :;n r f .ice boundary conditions t. a be prooaqat ed into 
the flow i iold ;o>nei, n'iiilt inq in better numerical behav- 
ior. in tiiis .;*• ud / , application of the boundary condition 
after t. he corrector resulted in instability, since a continu- 
al expansion was eventually attained, which dcove the body 
surface pressure to zero. 

The boundary conditions for the elliptic cone were de- 
rived in the physical Cartesian system tnd transformed to the 
computational polar system. This resulted in the application 
of the elliptic boundary conditions to a circular cone. 

Since the numerical method discretized the boly, the method 
could deduce the body geometry only through the boundary con- 
dition relations. 

Method of. Solution 

The coordinate system presented was selected primarily 
because it yielded a circular cone in tha computational sys- 
tem. ^ h«. significance of this, apart from the previously 
mention**.; inversion simplifications, was that boundary condi- 
tions could be tpplied to a constant coordinate surface 
( constant circular cone angle, 6 ). 
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The desired elliptic body shape could be acquired by 
choosing appropriate values for the semiraajor and senininor 
axes, a and b, and the characteristic length in the body axis 
direction c, which governed the slenderness of the configura- 
tion. After transforming the elliptic body to the computa- 
tional polar system, the number of angular increments around 
the conical body { N BETA ) were chosen, and the outermost 
shock position ( SHOCK ) was estimated. These choices de- 
fined the mesh increments, A(3 and Aft , used in the numerical 
integration, as follows: 


A(3 = 


130° 
N BETA 


(SHOCK - Cone Semiapex Angle) 
N DELTA 


(16) 


!* 


B. 


L 

£ 


Integration was performed from p * 0® to p = 180° in the 
positive sense of trigonometry ( see Figure 1 ) . As indicat- 
ed before, the 0 = 0 U , 180 1 * meridian is the plane of flow 
symmetry and therefore reflection across this plane was used. 

This reflection compensated for the backward and forward 
differencing requirements of the HacCormack predictor- 
coc rec toe seq uence . 

For the reflection boundary condition, the predictor and 
corrector steps were completed before the boundary condition 
was applied. This approach demonstrated good stability char- 

•1 
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.cteristics and ,1.1*.- satisfactory -Lti— “ 
lf noted. the entrap, boundary condition ... appUed after 

the predictor step. 


program execution was 
convergence criterion was me^ 


terminated automatically when 

This criterion was defined 


as: 


P 

,n-l 


( 17 ) 


„n*l i „ the pressure of the 

(here P n is the present pressure, 

a.- ..<* € is the convergence tolerance ( usu- 

ctevious iteration, and « is cne 

ally e = 0.001 ). Equation 17 was satisfied for ten consecu 
ti ,„ iterations before the esecution .as ter.inated. This 
gave rise to a conservative deviation of 0.1* P« iteration. 

The prog tan logic and u.derlyi.g theory .ere first 
tested h, investigating circular cones at rero and s.all 
indices, agreement vith other accepted solutions , 

bento, et al. <h) ) «ch that the approach .as dee.ed 

valid, and the elliptic con. prohle. ... then confronted. 
..suits of the couical coeput.tions are reported xn the nest 
section and the elliptic cone. nr. treated i. the sections 

followinq. 


THE CIRCULAR CONE 


In or lor to test the liability of the approach in the 
previous section, the reflection and entropy boundary condi- 
tions were applied to a 6 = 10° pointed, circular rone at 
- 5.0 and angles of attack of a = 0° and 5°. This circular 
configuration was selected for a number of reasons. First, 
in the past two decades there have been many circular con« 
studies, so that much comparison data ware available. .Sec- 
ond, the conical flow assumptions made concerning the flow 
field around th" elliptic cone were also valid far the circu- 
lar cone. Thus, the same approach could be used for both 
problems. Third, having a constant coordinate sjrface, the 
circular cone was mathematically simple to describe and an 
easy body to which boundary conditions could be applied. Fi- 
nally, some measure of merit was necessary to "judge the qual- 
ity of the reflection and entropy boundary conditions. Th*-"' 
circular cone, with its Dlethora of comparison material, sat- 
isfied these necessities. 

The initial concern w nr, lor the numerical location of 
the shock. Due to the finite difference mesh, the shock 
could, at. be st, be captured between two mesh points. Because 
of this, the step discontinuity in the d *rivativ»s of flow 
variables across the shock could not be represented exact iy 
by the numerical method. Therefore, there was so*° doubt is 


to whether th^-* shock position could be a 'CJ rat elv d !i pi'ted. 
However, as s °°:i from Figures 2,1,4 and 5, th» shock position 
is clearly defined and reasonably accurate. Figure 6 shows 
the numerical s.iock position overlying Babenko’s (4) compari- 
son data. The shock positions calculate! by e.aci bounlary 
condition were juite close and are represented by the same 
data point in the figure. As can easily be seen, the numeri- 
cal shock position representation is quite good. For the a= 
0° case, the angular shoe* position, from NAC.A 4-11.15 (14), 
is 6 -• 15.70, w |,ii t » t he numerical angultr shoe* position va 

calculated to b > 6 - 1>.o°. 

For the circular cone at zero anql? ot attack, th> en- 
tropy boundary condition invariably yielled a better repre- 
sentation of the pressure and density variation outward from 
the body. Figures 2 and .1 show that th? entropy boundary 
condition yielded pressure and density solutions that match 
those of Babenko, while the reflection bounlary condition 
yielded solutions that were 1.5^ to 2.0 % high. Although ob- 
viously not as u:c ur at •» as the entropy bounlary condition, 
especially at t so body surface, t.hv* reflactioa hour, lary con- 
dition should not be completely discounted. ">u-» to slow con 
vergence properties of the reflection boundary coniition, an 
a premium an computer time, the numerical intoqration was 
teruinated before convergence. However, as integration con- 
tinued, t he reflection bounlary condition solutions up- 



v , • i '*£*#&'?* ^*r*0i*. * 


proached tho.se obtained by Babenko. 

For the circular cone at angle of attack ( a = 5° ) , 
these same trends were witnessed. The conical body at inci- 
dence gave rise to pressure and density gradients around the 
body. This effect can be seen in Fiqur?s 4 an 1 5, where 
pressure and density are plotted for three angular positions 
( P - 0°, ‘*0° and 180° ) around the body. Tha windwarl side 
( p - 180° ), the region of greatest conpression, is aptly 
depicted by the outward pressure variations from the body to 
the crisp numerical shock. The leeward side of the body ( p 
= 0° ) , the region >f l«ast compression, shows a somewhat 
poorer pressure and shock representation. This is due to the 
eigenvalue structure and its previously noted effects on the 
numerical integration around the body. 

Of greatest concern were the density anomalies at the 
body surface, seen in Figures 3 and 5. This misbehavior re- 
sulted from the assumptions made regarding th? the flow near 
the body and the subsequent reflection boundary condition. 

To eliminate this problem, the entropy boundary condition wis 
applied after the numerical process had been brought to neat- 
convergence using the reflection boundary condition. This 
led to the most dramatic evidence of the beneficial traits of 
the entropy boundary condition. Figure 7 shows th*» efi<="~t of 

the entropy boundary condition. For the a - 0° case, th~ 

reflection density solution at the body was low oy 26^, but 


M 


.£. aw* 
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when the entropy boundary condition was used, the body value 
of density coincided with Babenko's rest Its. This same r.o- 
havior can be witnessed for the Ct = 5° case, although 
agreement at the body i? not exact. 

It was round that for coni iqurat io 1 s havir.i gradients of 
flow variables in the lateral direction ( the 0 -1 irect ion ), 
the numerical solution demonstrated a "lag" when the entropy 
boundary condition w.is applied. After only a f°w integra- 
tions ( 1 or 2 ) utilizing the entropy condition, solutions 
on the bod v vpia very close to the comparison data. However, 
as integration continued, further converjence was painfully 
slow. 

The introduction of the entropy boundary condition did 
not at feet the body pressure as much as it di! f he body den- 
sity, since the reflection boundary conlition yielded a nod 
representative pressure values. From Figure H , for the a = 

0° situation, the entropy boundary condition aqain brought 
the surface pressure solution into accordance with the com- 
parison pressure solution. ifor the a - = >° case, th=» charac- 
teristic lag was again encountered. 

These preliminary results led to th ■» conclusion that the 
entropy boundary condition better repreu»ntol th? actual flow 
field foi the circular cone. Therefore, iu» to the af>: •■'men- 
tioned similarities between the elliptic and circuit! cone 
problems, it was concluded thit the entropy bouniury cinli- 


*» '■**&&*. 
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t ion would yield the better representation of tha flow around 
the elliptic cone. Uith this Redaction, the plan of attack 

on the elliptic cone was formulated. 

The geometry and orientation of the confiscation would 
be defined, and numerical integration, using the reflection 
boundary condition, would be used *o acquire a nearly con- 
vorqed solution. This solution would then be compared to ex 
istin-j data and, it de^ed satisfactory, would be used as 
initial condition.* for rnrther numei ical integration, incor- 
porating the entropy boundary condition. This integration 
would :>“ continued until conver jence ( or aea « onveiqencr l 


was dt . 
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THE ELLIPTIC CONE 
Introduction 

A m.i )or problem in study ing the sun or sonic flow ibou t 
elliptic com':, is the me ,14 or supply of a ccessi bl » • j lli:)*ic 
cone conip.it i .or data. tot It analytic dat r and sufficient ex- 
peri mental information ,1 r 1 difficult to acquire. However, 
with the available liter at ire and consideration of the 
elliptic-circular cone similarities, a "thumbnail sketch" of 
the elliptic cone oroblen can be divined. 

Melnick (Id) indicates the presence of vertical singu- 
larities in t h-‘ semiminor axis Diane, at. the surface of an 
elliptic cone at incidence. According ta Melt irk, these vor- 
tical singularities, atisimi from the in vise il flow assump- 
tions, are the terminal points of th« crass flow streamline.. 
Helnick t hi-'n tone 1 ud es that "entropy, density and riditl ve- 
locity are Generally multivalued at the;^» point ;." floret ti 
(20) states that th« absence of viscous effects not only 
limits the validity of the solution, but also requires the 
angle of attack to be less than the conical semiapex angle, 
if a leeward shock is to be observed. Tn Sears (22) , F^rri 
calls attention to the possible lift-off of the leeward vor- 
tical singularity as angle of attack is increu«?l. sin.-*- th«» 
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angle of attack was kept small, the vortinl lift-off 


was riot 
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w i * ru-s scd . 

rot Mu' elliptic on? at zero inei i mi .:-», th» st a jni f ion 
at t ?a« line through the shock an: impinges on the *-l- 

liptic cone in the p - : J0° plane. Thus, the anti.ro elliptic 
surface, with the possible exception of the vortical singu- 
larities, is an i son tropic stream surfat'*. Since the loca- 
tion of the streamline impingement is known for the zero in- 
cidence case, tsuir* i i ca 1 calculation of the holy entropy is 
facilitated. liowevoi, for the anq] e-of- at t ack case, the 
stagnation streamline strikes the elliptic holy somewhere be- 
tween P '>0 ° a;M M0°, depending on the flow velocity and 
the angle oi at tack. Tie exact location of th ? impingement 
is pin-pointed when o the cross flow velocity is zero. 

To find the entropy on the body, the entropy of the 
stagnation streamline must be calculated Put sine? the nu- 
merical technique discretizes the flow field, the stagnat. ion 
streamline away from the body can only b? approximated by 
some interpolation method. Thia may leal to moleraf° errors 
in sui face entropy value:; and is a source of further study. 

Anothn otigii of probl.»ms is the area near t h •> - to 0 

plane. In this pen ion, t h»» large gradients of the flow viri- 
abl*-.; are dittioult to calculate precisely, duo to the dis- 
cretization or the t’.viw field. This behavior could, quite 
possibly, give rise to the characteristic M lag M noted for the 
circular cone solutions and accentuated in th» ollipti* con« 
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solution:!. A ?!;>re accurate description of tho flow field in 
the p * ( )i)° region roull be acquired by further clustering ot 
data points in this area. 

In o! ’er to facilitate the comparison with Chapkis' (7) 
result ;, the rollowing parameters were chosen for the ellip- 
tic cone at zero incidence: a/b = 2, c - 0 . 66 b 7 and = 5. 

For the ellipti: cone at angle of attack, comparison with 
Mar tel luce i ' s (17) theory required an elliptic cone with dif- 
ferent defining parameters. These parameters were defined as 

follows: a/b - 1.768, c = 1.0 and M = 3.09. 

00 

Zero Angle of Attack 

As in the circular cone study, firtt cons i l :-»i at ion was 
given to the shock position around the body. However, due to 
the lack of antlytic.il corroborative material, some amount of 
intuition must be used in "judging the numerical results. 
Chapkis (7) indicates that the shock should be closer to the 
body in the see ima jor axis plane than it is in the semiminor 
axis plane . This is due to the high pressure region near the 
elliptic "shoulder” ( i. p. , near the P - 90® plane ) and to 
the resulting pressure gradient iroa th'* .'haulier towarl the p 
0 ® anl 1 60 ° plane;. This observation is verified in 
Kiguie M. The shook positions calculated by *arh boundary 
condition were essentially coincident. Maximim pressure wa; 


noser ve ! 


90 ° 


.it st ign.it ion streamline, 1 k:U?1 on the p 

plane, which i>: floated the presence of a stvonq shock. As 
mention.'! before, the entropy value reqiired by the entropy 
boundary condition was calculated using flow variables ac- 
quired from this plane. 

The variation of surface pressure coefficient around the 
elliptic body is present el in Fiqure 10. The reflection 
boundary condition matches well Chapkis* (7) theoretical 
data, except in the range 80° < p < 100 °, which rouqhlv coin- 
cides with the region of greatest body curvature. Th° solu- 
tion obtained using the entropy boundary condition, also 
shown in figuio 10 , is far from the theoretical expectations. 
The figure indicates a more constant surface pressure coef- 
ficient for 1° <- 3 •" °0°, while the 70° < p < -90° region dem- 
onstrates a sharper surface pressure coefficient, gradient. A 
weak local expansion in this latter region woul .1 explain this 
behavior, although this is pure conjecture. (This will be 
discussed in greater detail in the follawinq section.) Due 
to the symmetry of the flow, the region 99 3 < p < 180° should 
match the 0° < 0 < region. However, Figure 10 shows an 

anomaly near P 1 1 ‘ » 0 f>r both boundary conditions. At this 
point, either ium.Mic.il "overshoot" or eigenvalue problems is 
suspected. Du.' to ' .it: k o f computer funling and t.im*, i more 
tio tough analysis <*o«) 1 d not be carried aut. 
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Angle of Attack 


For the angle-of-at tack case, the shock position around 
the body was again investigated. The results are presented 
in Figure 9. The point at which the shock most closely ap- 
proaches the body indicates the strong portion of the shock 
and the approximate location of the stagnation streamline. 

As the angle of attack was increased, tho stagnation stream- 
line moved in the direction of increasing 13. For a = 6°, 
the location of the stagnation point on the body was £ = 9B°. 
The uncertainty about the lateral stagnation point, and the 
extension of the stagnation streamline into the flow field, 
were major sources of difficulty in the application of the 
entropy boundary condition. Without a reasonably accurate 
approximation of the stagnation streamline position, entropy 
on the body surface could not be accurately calculated. 
Therefore, the boundary condition was lass reliable. 

Figure 11 shows the variation of surface pressure coef- 
ficient, calculated from each boundary condition, around the 
body, compared to »la rtel luce i • s (Ifl) theoretical data. Gen- 
erally speaking, agreement is good; however, the so° < 3 < 
riO° region is one notable exception. 

It seems intuitive that supersonic flow around a corner 
( such as the flow t ro* the stagnation ooint around the el- 



K. 


liptic should or ) would experience a sliiat, but not ice * 1> 1 e, 
expansion. No theoretical evidence of this was found in th“ 
literature avaiLaole; however, experiments conducted by 
Zaxkay and Visich ( 26 ) for a similar elliptic cone subjects 
to similar flow conditions, but at a = 150, s ^ em to Vf , L - ify 

this, conjecture. 

Figure 11 appears to indicate a locally expand'd flow 
for both boundary conditions in the previously noted b0° <p 
< 80 ° region, althouqh the entropy boundary condition yields 
a more pronounced behavior. In Figure 12 , the density around 
the body is plotted, showing the entropy and reflection 
boundary condition behavior. No comparison data was avail- 
able, but the behavior seems to lend credence to the “expan- 
sion hypothesis.” 

The validity of the underlying theory and assumptions, 
and the applicability of the entropy boundary condition for 
more closely approximating surface flow ch aracterist ics, have 
been demonstrated. flore thorough investigation of the numer- 
ical behavior of the entropy boundary condition and a noro 
concise method for obtaining the surface entropy are needed. 


28 


acknowledgements 


Th o to «pr« H. appcec iation to DC 

Dal e ».«». <« «. — - — " " 

. tul> and to th. I... «... a.t..»itT « — ^ »«•«= 

institute and Department of Aerospace Engine* for frnan 
cial support . The author also — to convert 
thanks to bis «if». Connie, for her patience th.oug 
period of his graduate education. Without her — 
and dauntless optimism, this graduate degree would h 

most difficult to obtain. 



23 




bibliography 


1 . Abbott, N » , Quarterly '^nort , Aerot.h?rn, Proiect 61«7, 

Aug . 1971. 

2. An dor. son , 0. A. and Vogel, J. ri. , "Numerical Solution.'; 

!• 1 i>wt if Ids IV h ind Rectangular Wings," Engineer ing i)o- 
searoh Institute, Iowa State University, Am«s, Iowa, duly 


1. 


a. 


6. 


An lot son D. A. and Voqol, j. n., 
of Flowf i -'Ids, " o uarter 1 y Report , 
I n s t i t u t *■ , low » Stat<? University. 

1971. 


M N m nr i o * 1 Cal o u 1 i t ions 
En jineerin-i Research 
Am^s, Iowa, July- Aug. 


Babenko, V. I., Voskresenskiy, G. P. , Lyubimov, A. N. and 
Pusanov, V. v., "Three Dimensional Flow of Ileal Gas Past 
Smooth Bodies, " NASA TT F-dRO, Apr. 1966. 

Briggs, t< . R . , "The Numerical Calculation of Flow Past 
Conical HodUs Supporting Elliptic Conical Shock Wav^s at 
Finite Angles of Incidence," NASA TN 9-340, Nov. I960. 


6. Biook, j. w., "The Calculation of S>nlinear Supersonic 
Conical Flows by the Method of Integral Relations," ? light 
Dynamic ; Research and Technology Division Report n 4 - 7 . 

July 1964. ' 

7. Chapkio, F . L., "Hypersonic Flow ov or an Elliptic Cone- 
Theory and Experiment," Journal of th* Aerospace Sciences 
Vol. 29, No. 11, Aua. 1953, rp. 844-964. 

B. Ferri, A., "Supersonic Flow Around Circular Con»r, at 
Angles of Attack," NACA ^-2236, Nov. 1959. 

9. .erri. A., "The Linearized Characteristics Nathod and 
its Aoplicition to Practical Nonlinear Supersonic Prob- 
lems," NACA TN- 2515, Oct. 1951. 

10. Ferri, A., Ness, n. and Kaplita T. T. , "Supersonic Flow 
over Conical Fodies without Axial Symmetry," Journal of 
the Aerospace Sciences, Vol. 20, No. 6, Aug. 1951, on. 
r )h r ,7 1 . 

11. :!oid, P. A., "An Approximate solution for Axially ;ym- 

:v*M ic Flow over a core with an Attache? Shock Wave," vn-» 
TN- 34M r > # oct. 1-56. ' 


■j 

A 


30 


12. 


1 A, 


1 U. 


is. 


if.. 


17. 


Ifl. 


19. 


20 . 


2 1 


,f t ',o PUw c i^l r ' ^ oc 
r , ! ..Nu*oricdl Solut "°”Lr- National Ae ronau- 

L plow About 

f ind c f>ldl^kl # £ - on •: Journ ll ot n 

Kihanv'# A. ) n<1 pniptio Cross section, * 

po-i.o* ; »o. 8. Wj ’. n 

,eL0Si..c- sn«™ ,J - vo '- 

. . i) if f o r a nCt- J 

...K, P-. " k f l ^Uto,°o£ S C^icU «o^jjobl.«.* 

Oct. 1989. . . 0 j- Moffett 

-uur.: sssr.s «-» - 

'‘°" H , , P "Tb 3 Numeriv ii so 

Jt‘r- -,1£. -•» 

V-7 69# r iM/ in , lYl^ rVO ' 

. t v, ••The Effect of Vl ^’ 4-3SU, 1999, P»* 

,v. ;:Cor«aek , ' •• . „ kI y Paper No. 

locily Ufict Ciatai m j, 

1_7 ‘ of th- lin^riw' Charac- 

sr'«! pp - 667 ' . , 

<n conical 

,elnlck. F. "I OC ;o!! l S l So- '«<• 1,67 ' P! ” 

Flow,” AT A A Journal, Vol. - 

631-637. . _ « mv« "low," 

■skVa. Essarass&’a* >./- - *•- 

noretti, <••» l "* KlhK journal, vol. # 

,»t an An<|l*’ ot Attack, 

vjk7 up. 

1?b7 ' Fl ^ , Ko of Ch » i act.er- 

«.K»CF. .. * V-V rtf-'AuZ*' 

nia. 


22 . 






2 3. Spars, W. P., ed., General Theory o£. Hi ail SEged Aerody- 
namic;;, 1st ed. # Vol. 6, Princeton University Press, 

Pr inept. on, New Jersey, 1965, pp. 721-746. 

24. South, J. C. and Klunkor, E. 3., "M»thols for Calcu- 
lating Nonlinear Conical Flows,” SP-226, Oct. 1969, Ames 
Research Center, Moffett r iold, Cali r ornia. 

25. Van Dyke, m . p., "The Elliptic Cone as a Model for Non- 
linear Supersonic Flow Theory," Journil of Fliud Mechan- 
ics, Vol. 1, No. 1, 1956, pp. 1-15. 

2 6. Ward, G . N . , "Supersonic Flows Past Sl?nd j r Pointel 

Bodies," Quarterly Journal of Mechanics and Applied Mathe- 
matics, Vol. 2, 1949, pp. 74-97. 

27. Zakkay, V. and Visioh, M. , Jr., "Experimental Pressure 
Distributions on Conical Elliptical Bodies at M <*, = 3.09 
and 6.0," Polyt echnic Institute of Brooklyn, PIBAL Report 
No. 467, Mar . 1959. 




ndiE co 






lar Position Outward from the Body, 






Angular Position Outward 











*HL, 



Xfcfc 



Figure 3 Lateral pressure 












Chapkis; 



Lateral Body Angle, 




V»-t . ’ x :' -T* ’>/♦>. 



Lateral pressure distribution 







44 


APPENDIX A: COORDINATE SYSTEM AND CONSERVATIVE EQUATIONS 


Coordinate System 

The coordinate system chosen, Fiqure 1, required two 
transformations to acquire the final computational system. 
The first transformation was from the Cartesian system, 
(x,y,z), to the computational Cartesian system, (x',y',z')» 
via the transformation: 


x ' = x 

y* * ay 

z* - bz , 


( A - 1 ) 


where a and b are the semimajor and semiminor axes of the el- 
liptic cross section, respectively. 



Transfoxmation A-l 


Fiqure A- 1 Coordinate system transformation 


The effect of the transformation proviiei by Equations 
A-1 is the def ormation of the elliptic cross sections in the 
Cartesian plane into circular cross sections in the computa- 
tional plane, see Figure A - 1 . 

The general equation for a right elliptic cone is: 

<V + (V ^ (*/ • (A '' 

b a c 

Referring to Figure A-1 and Equations A-1, it can be seen 
that Equation A-2 can be written as: 

< y ') 2 - ♦ (z') 2 - . < R ~ 

C 

The final transformation takes the computational 'arte- 
sian system into the computational polar system, (?,6,0), b 
utilizing the geometric relationships in Figure 1: 
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Conservative Equations 

The MacCormack second order preferential predictor- 
corrector requires that the equations of fluid flow. Equa- 
tions D- 6 , D-7, D-8, D-9 an 1 D-10, be written in 
conservation- lav form. These fluid flow equations can be 
written in matrix fora and represented by: 

E, ♦ F y ♦ G z = 0 => E x . ♦ aF y , ♦ b.l^ = 3 . (A-7) 

The desired fora is: 


E r ♦ aF. ♦ bG ♦ H = 0 

5 6 0 


(A-*) 


Froa the chain rule and algebraic manipulations. Equation A - 8 
becomes: 


Ep - ( n cos 6 sin 6 ) ♦ a cos? 6 cosp ) F^ 

♦ b( i (cos* 6 sinp) G ♦ i (cospcot 6 ) G J 


( i cots siip ) F^] 

0 . (A-7) 


With S # 0, i.e. , away from the conical apex, the conserva- 
tive fora is obtained as: 


fE Zip* f-Vsin 2® ♦ (aF cos9 ♦ bG sin0) cos*6] 

•' 6 

♦ f(-aF sinp ♦ bG cos9) ~ 3 t 6 ] (A- 10) 

9 

♦ f (co s 26 - 1) E ♦ (aF cos 9 ♦ bG sin 9 )(cot 6 ♦ sin 26)) 


0 . 
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Pro* Equation A-10, the eatrix repcesenataion of the fluid 
flow equations in the conputat ional polar systew becoaes: 


§( pu } 

■ 

5(P* PU2 ) 

A 

5 ( puw) 

T 

f(puw) 

F L 


{ pu ) sin 2 6 
-^(P > pu 2 ) sin 2 6 ♦ 
- \ ( pu tr) sin 2 6* 
-k ( puw) sin 2 6 ♦ 


♦ ( a(pv)cosp ♦ b(pw) sinp)co^6* 
(a(puv)cosp ♦ b( puw) si np ) cos 2 6 
(a (P* pv 2 ) cos p ♦ b ( pvw) sinp) co^6 
(a(pvw)cosp* bfp+pw 2 ) sinp) co^e 


+ 


(-a(pv)sinp ♦ b ( pw) cosp) cot 6 
(-a(puv)sinp ♦ b ( puw) cosp) cot 6 
(-a(P* pv 2 ) sinp ♦ b ( pvw) cosp) cot 6 
(-a ( pvw) sinp ♦ b (P+ pw 2 ) cosp ) cot 6 


3 


(A- 11) 


+ 


-2(Pu)sin 2 5' •> (a(pv)co.sp * b ( pw) sinp ) (cot 6 * sin26) 

-2 (P* Pu 2 ) sin 2 6 ♦ (a(puvjcosp ♦ b(puw)sinP) (cote ♦ sin26) 
-2(puv)sin 2 6 ♦ (a (P* pu 2 ) cos P ♦ b ( pvw) sinP) (cot 6 ♦ sin26) 
_ -2 ( Puw) sin 2 e ♦ (a(pvw)cosp ♦ b (P* pw 2 ) s inp ) (cot 6 ♦ sin26) 


-0 


which has the general tors of Equation A-8. 




7 




*. ■- »v. 


\ £ 
►. * 
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APPENDIX B: THE REFLECTION BOUNDARY CONDITION 


The equation of the elliptic cone in Cartesian coordi- 
nates can be given by: 


2 2 2 n 

♦ / z \ - t - 0 • 


P(x,y#z> * (i) ♦ (i> - (JL) 

b a c 


(B-1) 


Due to the conical flow nature and the fact that the body way 
be replaced by a strean surface, the exact boundary condition 
is: 


q* VF = 0, on the body, (B“2) 

so that the velocity relationship that nust be satisfied on 
the elliptic body surface is: 

v m - JL * ( Jl) 2 (B~3) 

u y c y u a 



An alternate way of expressing Equation B-2 nust be found, so 
that the numerical aethod can be used. Froa the geometry of 
the problea, see Figure 1, it can be seen that: 

* <? « V (B - U ' 




Fro» calculus: 



iq N i 


= i , ~u v tan6cos|3 . v tanosmp n 

3 ir* + — -h 1 *— “ + ,! h ' 


(B- 


ab 


a 2 b 


where 




(B- 


The reflection boundary condition chosen assumes: 


iq„i 4 = (P ~ 

and 

|q | ♦ * I q |~ , (B- 

14 tan ' 1 M tan 1 

whe re the (♦•) sign indicates the values at the first mesh 
layer above the body ( superlayer ) , and the (-) sign indi 
cates values at the first mesh layer below the body < sub- 
layer ). This component reflection can be visualized as: 
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tan , tan 


yy 7 ■■ r j 1 / r r rr y 

• i 





Super layer (+) 
Body Surface 
Sublayer (-) 


Figure B-1: Component Reflection. 

Equations B-10 and B- 1 1 actually yield an overdeterm ined 
system. Since chere are only three unknowns (u , v and w ) , 
only three reflations are needed. Let: 



q tan “ rj-j^tan *, # * q tan* 2 ' * q tan* 3 * 

t 

(B- 15) 

which 

has form identical to Equation B-10. 

The 

system of 

three 

equations necessary to solve for the 

t hre<= 

» unknowns is: 


‘V* = -%•' 


(B- 16) 


'W* = 


(B- 17) 


*^tan ” ^tan *3 


(B- 18) 


Thus from Equations B-10 and B-11 the following is acquired 
for the boundary condition relation: 



where the (♦) superscript indicates that the variables are 
evaluated at the superlayer and the (-) superscript indicated 

evaluation at the sublayer. 

Equations B-19, B-20 and 3-21 constitute a system of 
three equations and three unknowns which, when written in 
aatrix fora, becoae: 



* a 


■t- 

*; 

i 

I 


I ' 


• 


i 

1 
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- 1 q„ I A 

S 


* ^ tan ^ ) 


L 


"'tan 1 , 


{ B- 22) 


eraser's rule can now be applied to Equation B-22 to deter- 
■ine the desired values of u~, v-, and w~. 

Pressure and density are assuaed to be even functions 
across the body boundary so that their normal derivative can 
be set equal to zero. This yields the pressure and density 
boundary conditions: 


P* = P~ 


(B- 23) 


and 


P* * P 


( 8 - 214 ) 


r -~ 






T 



X • 



> ^5 

* 


From Figure C-1, it is seen that: 


— * A 

-4*_ ♦ i ' cos e = sin(AG) (C-1) 

I c. | N 


so that. 


A0 




(C-2) 


where *A0 denotes an expansion and -A 9 denotes a compres- 
sion. Carrying out the operation indicated in Equation C-2 
yields: 


-fiii 


A 0 = 


- I 

sin [ £ 


v, Y 


iq, i * 


(C-3) 


where d is defined by Equation B-9 and i N is as in Equation 
B-6 . 

Fro« NACA B-1135 (15), the static pressure after turning 
the flow parallel to the body is determined by: 

Pa =P, [ 1 - y«i-(A9) + -Hi 2 ( v* 1 1 H* - 4(H, 2 - 1) (A9> 2 . 0 (A9 ) \ 

W- D* 4 (H, 1 - 1) * J (C-4 ) 


where M, is found from: 



56 


M< = 


|q. I 

a. 


iq. I 


* [ 


(u, 2 + 


v, 2 + 




-*v 


(C-5) 


Note: P, , p, , 1,, | q . |, and a, are all nonditnensional vari- 

ables as defined in Appendix E, unless otherwise specified. 


From the equation of state, p 2 = P 2 (S,P 2 ) t that 


_ p 2 _ constant 

P? ’ pT 


(C-6) 


which yields the relation for the density: 


P 2 



1 

y~ 


(C-7) 


where P t and Pj are obtained from the plane containing the 
crossflow stagnation point. For zero angle of attack, the 
crossflow stagnation point coincides with the p = 90° plane 
( see Figure 1 ) . 

From the energy equation, total enthalpy is found as: 


h 


T * 





+ 


iJj 

2 


2 


(C-fl) 


which yields the velocity modulus, |g 2 | : 


iq* l 




(C — 9 ) 


Note: In Equation C-8, P , F, and |^| are in dimensional 
form. 


Thus, the magnitude of the velocity vector tangent to the 
surface is found. In order to determine q 2 , a unit vector in 
the direction of q, is defined as: 

iq £ Uj'Vj 1*2 , w| . cc- 10) 

2 iq 2 i fgTi 


rrom Figure C-1, 

3 2 * 3 tan = 3 ; = ?, - ( 3 . • i N )i N • (C-11) 

Therefore, q® is in the direction of q^but it is not neces- 
sarily true that |qj | = |q 2 |, due to truncation error. The 

A 

unit vector, i a , can now be formed as: 
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Carrying out the operations indicated in Equation C-11 
yields: 



1 f Jh-JL 2 _ Yt *y _ w. xz . 
d l v c* b* c l a J c*' , 


' * ' + V * N*>>. 


[ 


w. 


J / -IUJLZ 
d* ( a 2 c 2 


Yi yz , 
a^b 2 + 



(C- 13) 


where d is es in Equation B-9. Finally, 


3 » = 


iq 2 1 


* 

-12. * Iqa I x q 
If; I q 


(C-14) 


can be calculated. 

Since numerical computation is performed in the computa- 
tional polar systea, the appropriate equations above most be 
transformed using Equations A-1, a-4, a-5 and A-S. This 
yields: 
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and since q® = (u®,v°,w°): 


u * = U ’ “ tanfisine . (c-16) 

M m c 4 ab 2 c l a l bc l 1 


v 2 g V,- (Jw l li- ta y 6yO§& -tt- taD ^CP5 2 P - * tan a 6 ^Q^6sinfl ] (C-17) 


"* " w « "('ll 2 f tt i- tabSsjp fl _ Yi__toLn 2 fi cospsinp _ i, tan 2 fi sinp i (C-1R) 
d a* be 2 a’b> ^b~2 UJ 


where d is froa Equation B-12. Finally, 


*1 



Hrfr ( ^r < 1 


Pa 


)] 


(019) 


and likewise. 


v a = _lL Ha I and w, 

\W\ 



(C-20 ,21) 


The crossflow stagnation points are located at the 
points where q c = 0. The derivation for the radial and 
crossflow velocity components ace now presented. 

Froa Figure 1, it can be seen that: 





D = (x* ♦ y 2 ♦ Z 1 )* . (C- 25) 

Proa Figure 1, it can be seen that: 


and 


i«c i - <r • 

A 

‘c 

(C- 26 ) 


A 


*«R ' = 1 • 

l R 

(C-27) 


However, since numerical computation is again performed on 
the computational polar system, the transformation equations. 
Equations A-1, *-4, k - 5 and h-6, must be appliel to Equations 
an ^ C-27. Carrying out the indicated operations, Equa- 
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APPENDIX D: EIGENVALUES 

The matrix form of the governing partial differential 
equations of fluid flow in the Cartesian system is written 
as: 


Au ♦ Bu ♦ Cu =0 . (D-1) 

§ « e 

However, to acquire a form amenable to the HacCormaok second 
order predictor-corrector technique in the computational 
polar system. Equation D-1 must be rewritten into the form: 

u ? ♦ A” Bu ft ♦ A"Cu p = 0 , ( D — 2 ) 

where A, A, B, B, C and c are nxn matrices, n being the order 

•« T 

of the system, and u = (u,v,w,P,p). 

To derive Equation D-2, transformations must be made 
from the Cartesian system to the computational polar system 
using the chain rule and the transformation Equations A-1, 
A-4, A-5 and A-6. In this manner. Equation D-2 becomes: 

u, ♦ [ Ar ♦ a A** B 6y ♦ bA^CA^tl. ♦ IaA"*B^ v . ♦ bA^Cp^JII = 0 , 

S Op 

(D-3) 


where 
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A = A 

B = A 6^ ♦ aB6y' + bCA* 

and 

C = aBpy/ ♦ bc Pr' 


(D- 4) 


fc-iR and A -1 ? «ust be found. The pac- 
Now the lattices A *B and * 

tial differential equations of fluid flow: 


X-aoientun : 

(P ♦ pu* \ 

♦ ( pUV)y ♦ ( PUW) I 

= 0 

(D-S) 

Y-Ho#entu»: 

(piiv), ♦ 

(P ♦ pV 2 )y * (P V- >I 

S 0 

( 0-6 ) 

Z-rtonentus: 

( 0 uw) x ♦ 

(pVWjy ♦ (P ♦ P W *>Z 

= 0 

(D-7) 

Energy: 


q . (TP - cVp ) 

= 0 

(D-8) 

Continuity: 


T * ( p3> 

= 0 , 

CD— 9 ) 


where free Kutler (14), c* 

B and C. The desired quanti 


- (d_L)| , yield the satrices 

dp S 

ties are then obtained as: 


r. 



c 2 -u 2 


c« -u* 

CM 

2 

1 

M 

n 


pw 

0 

;p u 

- w 

w 

c 2 -u* 


C z -u 2 

u ( c 1 - u 1 ) 

u 


utilizing Equations A-1, A-4, A-S and A-6, Equations D-4 can 
be rewritten as: 

A = A 

3 - - J_ cos 6 si n 6 ♦ a A -1 B cos*6cos0 ♦ b k* C cos a 6sinp 

? ? ? CD-13 

C = -a K 1 B cotfisinp ♦ b A* C cotficosp 

5 5 

The eigenvalues, and eventually the linear stability 
bounds, are found by sol /ing the characteristic equations: 

det ( 0 - ct I ) = 0 (D-14 

and 


det ( C - XI ) 


0 


(D- 15 



and obtaining the characteristic roots. The resulting eigen 
values are: 


ct = -cos6sin6 ♦ cos 1 6 (avcosP ♦ bvsinP) 
u t » 3 | u | 

ct # 5 = [-((avcosp ♦ bvsirp) ucosfi «• sinft] 
’ «• (C*-tt k ) 


+ { [ (avcosp *■ bvsinp) ucosft ♦ sin6 ] 2 
“ ( <c 2 -u 2 ) 


sin*6 -2(avcosp+ bwsioS ) u cos 6sin6 

( c* - u * ) 


♦ (avcosp ♦ bwsif>9) 2 cos 2 6 

(c r -u r ) 

L 

-(a* cos 2 p ♦ h 2 sin 2 B ) 0 *00 5*6 7 | cos6 

(c^«*j 3 J 5 

and 


X = (bwcosp - avsinB ) cot 6 

»» it 3 u $ 


X ~ (avsinp- bwcosS) u 

M l (C 2 “U 2 ' 


+ { (avsins - bvcose) 2 u 
“ 1 <c*-u 2 


) 


♦ (avsin3 - bwcosflf 1 


-(a f sin 2 9 ♦ b > cos , 8) s 2 7* 

(c* -u 2 ) J 


cot ft 

t 


(D- 16 


(D- 17 


(D- 18 


( o - iq 



t>7 

appendix E: nondimensionalization and initial conditions 

Nondi«ensionali 2 ation 


The dinensional flow variables P, p, 
nond i Bensiona li zed to P, p , u , v and w as 


u, v and w are 
follows: 


P = 


P 

VP.t 


(E-1) 



(E-2) 


u#v,w 


[ lls: Ti* 


] 


(E- 3 f 4 r 5) 


Initial Conditions 


Fro» 

NACA 

R-1135 (15) ; 







-v 





_P_ = 
P ocT 

i 1 * 

(V- 1 ) .1 * => 

2 1 

P = 

It 1 

Y 

~V 

♦ (v- 1) wl 1 ^ 
2 J 

(E-6) 

P_ s 

Porfj 

t ' * 

"1 

2 J 

P * 

l 1 

-l 

(E-7) 

-i. = 

q« 

i ’ ♦ 

=> 

2 1 

q * 


♦ UiJ) % 

2 1 

(E-fl) 


since 






(K-9) 
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APPENDIX F: INVERSION 


Fron Appendix A, Equation A— 1 1 # define: 


E, = 
E a = 
E, = 
E* = 


5(p ♦ pu a ) 
Spuv 

IPUW . 


(F-1) 

(F-2) 

(F-3) 

(F-4) 


then the flow variables can be found as: 


E, 

E, 


(F-S) 


v * F.f 


(F-6) 


Ei 


U = E, - f E> 


4e; 


( Y-1 ) lr I . (V-1) (V* ♦ v*K7 

2Y 1 1 Y 2y J 5 

— in? rn 


(F-7) 


p « EL 


Ej_ 


(F-8) 


P = E a - Pu J 

> 


(F~9) 


The (♦) sign in Equation F- 7 is used because the flow is 
supersonic. 
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NUMERICAL SOLUTIONS OF THE ADVECTION EQUATION 
USING CONSERVATIVE FORM 

Dale Anderson and Behrooz Fattahi 
Department of Aerospace Engineering and 
Engineering Research Institute 
Iowa State University, Ames, lova 50010 

Abstract 

Second- and third-order finite difference methods recently applied 
to problems in high speed fluid flow are applied to the model advection 
equation cast in ccnservative form. The differencing methods considered 
use forward time differencing with preferential space differences em- 
ployed in the predictor-corrector sequences. The free parameter required 
for stability in the third-order method is adjusted to cause the solution 
to be either minimum dispersive or minimum dissipative in nature* Results 
indicate that the third-order method using minimum dissipation is the 
most accurate method tested. Computer time requirements are approximately 
twice those needed for second-order techniques. 


1. Introduction 


Numerical integration of the equations describing atmospheric 
notion has become commonplace since the advent of the high-speed digital 
computer. Techniques used in performing this integration have varied 
widely in order of accuracy, programming ease and computation time 
required. This paper deals with application of recently developed 
second- and third-order techniques to the advection process using the 
differential equations describing advection in conservative form. 

The accuracy with which a finite difference method models advec- 
tion within the atmosphere is difficult to assess since the governing 
equations are extremely nonlinear. Molenkamp (1968) has suggested that 
representative numerical methods be applied to a simple problem whose 
exact solution is known in order to compare the results obtained through 
the use of each technique. The results of his investigation show that 
higher order numerical methods more accurately model the advection pro- 
cess. Molenkamp concludes that of the methods tested, the Roberts and 
Weiss (1966) technique most nearly approximates the exact solution. 

This method is fourth-order in space and second-order in time. The dis- 
advantage lies in the fact that it requires at least one more order of 
magnitude in computing time and considerably more storage space than, 
for example, a Lax-Wendroff second-order scheme. 

Crowley (1968) has performed a similar analysis including a compari- 
son of the accuracy of solutions obtained using advective form and conser 
vative form of the describing equations. Tnese results were obtained 
for longer integration times but generally agree with those of Molenkamp. 
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In addition, Crowley concludes that the conservative form of the governing 
equation is preferred over the advective form. 

The test problem which Crowley and Molenkamp have used assumes that 
advection of a quantity A is governed by the differential equation 


dA dA dA 

dt " " U dx " V dy 


( 1 ) 


where u and v are the velocity components in the x and y directions 
and t is the time. This equation can be put in conservative form if 
the divergence of the velocity vanishes as 


du . dv 
dx dy 


0 


( 2 ) 


The nonhomogeneous terms vanish and the conservative form of (1) becomes 


The conservative form (3) and the advective form (1) describe the same 
physical phenomena (assuming Div V * 0) and should provide the same 
solution for a given set of initial data. Since numerical solutions 
are obtained on a finite grid, the solutions to (1) and (3) are different 
except in the limiting case as At, Ax and Ay approach zero. 

A particularly simple model of the advection process can be studied 
when the velocity field is prescribed and due to a rigid body rotation 
of angular velocity 0 about some axis. In this case, the advective 
form may be written in polar coordinates ar 


£A 

dt 


u de 


o 


( 4 ) 




u 
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where 0 is the angular coordinate. This is just the simple one-dimen- 
sicnal wave equation which has the solution 

A(r, 6, t) *= A Q (r, 0 - fit) (5) 

where r is the polar radius measured from the axis of rotation and A Q is 
the initial distribution of A. The initial distribution is rotated 
through angle fit. This is the correct solution to either the advective 
form or the conservative form of the governing equation. 

Application of numerical differencing methods to this simple problem 
provides a rough comparison of the usefulness of the methods in modeling 
the advection process. Unfortunately, the velocity field usually depends 
upon t and is coupled complexly to the advected quantity A in the real 
problem. However, by solving the numerical problem in a rectangular 
coordinate system, some of the nonlinear character of the original 
problem is retained. 

In the following sections, several numerical techniques are reviewed 
which have been applied to high speed gas dynamic and aerodynamic problems 
with great success. Phase error and amplitude plots are constructed and 
the methods are applied to the conservative form of the advection equation. 
The results obtained are compared with those of Molenkamp and Crowley. 

2. Differencing Methods 

The numerical techniques used in this study have been developed for 
hyperbolic systems of the form: 



& 

dx 


0 
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where E and F may be N component vectors and x may be of arbitrary dimen- 
sion. The techniques are based upon the Runge-Kutta method and are 
explicit in advancing the solution forward in time. 

McCormack's Method MacCormack (1969) developed a noncentered second 

order method which has been used extensively in solving gas dynamic 
problems involving strong shock waves. This technique consists of a 
two-step predictorrcorrector sequence of the form 


E*j +l - E" - g (1 - eA)AF* 


*r l ■ H e " + * r l ] + 


Ml+l 


(7) 


( 8 ) 


where the tilde represents values at the intermediate step and 


* F J * Vi - F J 
- F J - Vi • 


(9) 


The value of e can be either 0 or 1. If e * 0, the predictor uses a 
forward difference and the corrector uses a backward difference on F^. 
If e » 1, the differencing in the predictor and corrector are reversed. 

This method was derived using staggered differencing with the con- 
straint that both left and right moving waves dc properly computed. 

In addition, the method satisfies the perfect shift condition. If the 
simple one- dimensional wave equation is considered 


3u s$u . 0 

at a* 


( 10 ) 
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a solution may be written 


u(x, t) = u(x • cAt , t) 


If v = , then it follows that 

Ax 


u(x, t) = u(x - I/Ax, t) 


(HA) 


For a finite difference mesh the perfect shift condition 


n+1 n 

U j ' U J±m 


1 / - + m , m — 1, 2, 


is satisfied where v is the Courant number. 

MacCormack's method car be viewed as a version of the Lax-Wandroff 
method. The use of noncentered differences is advantageous since only 
values on the mesh points need be computed. This eliminates difficul- 
ties encountered in storage shifting and results in a considerable 
reduction in computer time. Kutler, Lomax and Warming have recently 
extended this method so the predictor can be evaluated at any fractional 
time step. Only values at (n + l)At were used in this study. 

MacCormack's method is readily extended to two-space dimensions and 
when applied to (3) becomes: 

~ a T] - A ”,j - “ [ <i - c x> <uA, w,j • u - 2£ * )(ua >;.j 




U - e y X »*) li)+1 


U - 2e y )(vA)J ^ - e y (vA) 


M-l] 
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A 1J ■ 2 [ A ",J *K'i\ ■ 2 IS + (1 - 


. , i\/~i\ n+ i I . at T .n+i 

+ <e x - IXuA).^ jJ + s [€ y ®A) tj 


J+l 


k n+l 


+ (1 - 2e F )<DA). ' + (e - l)CA) 

y ^ > j ■* * y 




(13A) 


Since the e parameter can change in both the x and y direction, a ques- 
tion arises in selecting the best values of the e's to use. MacCormack 

(1971) suggests that the values of e v and e„ be cyclically permuted to 

x y 

obtain results which are spatially unbiased. The allowable values of 
the e's are: 


€ 


e 


e 


e 



(14) 


Only values e x - e y « 0 were used in this study. That is, only forward 
prediction-backward correction was investigated. 

The stability of this method was investigated in detail by MacCormack 
and his results show that stability is assured if 


M £ 1 (15) 

where v it the Courant number given by 

- AT 
I ;*ffr 

Ax 

and a represents the eigenvalues of th<. Jacobian matrix dF/dE which 
arises in the stability analysis of (6). The stability bound is derived 
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i 


by noting that the absolute vaLue of the eigenvalues of the amplifica- 
tion matrix must be less than one where the eigenvalues cf the matrix 
are given by 


- 1 - il/ sin 0 - V (1 - cos 0) 


(16) 


In the application of this method to the one-dimensional wave Eq. (10), 
the amplification factor of the scalar difference operator is given by 
(16) where the Courant number is written as cAT /Ax. 

It is instructive to consider the phase angle and amplitude of the 
one-dimensional amplification factor 

2 

% » 1 - iv sin 0 - i* (1 - cos 9) . 


The magnitude is given by 


|§| - [l - 1/ 2 (1 - l/ 2 )d - cos 0) 2 ] 


1/2 


(17) 


and the phase is given by 

A . fan " 1 Im( ^ 

0 tan Re(5) 


tan 


-1 


- i/ 8 in 9 


(18) 


l - 1/(1 - cos 9_ 


As noted by Crowley (1968), the relative error in amplitude is given 
by 1 - and the relative phase error is 1 - 0/Ov when the solution 
is advanced in time. Figure 1 presents contour plots of constant 
for values of v and the reduced frequency 9 lit. The contours shown in 
Pig, 2 are lines of constant - 0/8*/* These plots show that both phase 
and amplitude errors increase with 0/ tr for a given value of t/. The 
amplitude error is small at low values of |/. The amplitude error is 
••all again at large values of p< The phase error does Just the 


i 
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pposite. The amplification factor and the plots of ampl ^ 

lngU ate Che same as »ere preaented for the Lar-Wendrof - 

8 . for these two methods is Identical. 

Irowley. The amplification factor for 

. Me thod Rusanov (1,70) and Berstein and hirin (1,70) 

. . thltd . order scheme Based on application of the 
t.neously developeo a third t€clml , u . has 

Runge-Kutt. method. Due to It. recent deve ° P "^^ ^ McCormack's 

net been a. extensively applied to fluid flow P* (Wn) 

thod Taylor. Hdefv and Maason (1,72) and Anderson an 

meth d * y Kn w h Rusanov ' s and MacCormack s 

J ms.ratlve solutions using both Rusanov 
have presented comparative 

, a a mat the third-order method possesses 
methods. They have concluded that . £leld 

greater flexibility in problems where larger variations in 

r It should he noted that the results of those inves- 
properties oocu . relatively short compu- 

tlgations were for shocked flow, and represent 

8 U finu field dependent variables. 

ration time, with large change, in « (6) 

-hod aoolied to the hyperbolic Eq. (b) 

The third-order Rusanov method appl 

is 

E 


(1) 

j+1/2 

-H 

n . E n \ _ 

E J+l + j) 

( F " + i ‘ F 

b (2) 

B j 


3 Ax V 

F (l) . 

F j+1/2 


_n+l 

E j 

- E ” 

3 AT | 
" 8 Ax 1 

- 2F (1) 
L 2 J +2 

t - ( ;i - 



3 AT 
8 Ax 

[f (2) - 
L F J+1 





n - 4E” , + bE*J 
J+2 J+1 i 


( 19 ) 
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This third-order method was developed using centered differences. The 
perfect shift condition is also satisfied for a Courant number of one 
and the proper choice of the to parameter. The omega term in the third 
step is a fourth-order difference required for stability. Since it is 
fourth-order, the third-order accuracy of the method is unaffected by 
its addition. 

The stability analysis of this technique has been carried out 
in detail by Burstein and Mirin (1970). Their results uhow that the 
stability of the system is assured if 

M <1 

and 

4v> 2 - v* < (0 < 3 (20) 

The application of the Rusanov method to the linear Eq. (10) produces 
the amplification factor of the scalar difference operator which is 

% - 1 - sin 2 9 - ^ sin** j 

- iu sin 0 £l + ^ (1 - V 2 ) sin 2 (21) 

Contour plots of phase and amplitude are presented in Figs. 3 and 4. 

An additional difficulty arises in presenting this data since the factor 
to appears as a parameter free to take on values within the stability 
bounds given in (20). The only possible contours that can be constructed 
and interpreted are those for which to * 3. The lower bound on to pre- 
dicted by linear thoery depends upon the Courant number. As a result, 
the curves for any value of to less than 3 would violate the stability 
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bound over some range of |„|. Figure 3 presents the constant magnitude 
contours for |5|. In contraet to the second-order method, the largest 
damping occurs at midrange values of wave number while at both high and 
low wave numbers relatively little damping occurs. In comparing the 
second and third-order amplitude contours. It appears that little is 
gained by using the higher order. The phase error contour, show an 

Improvement . 

These results are better understood if the modified equation pro- 
duced by applying each method to the wave Eq. (10) U examined. Applica- 
tion of MacCormack'a method to the linear wave equation results 

modified Eq. 


du cdu 
dt dx 


- i cAx 2 C 1 - l/ 2 ) ^3 + ••• 
6 dx 


( 22 ) 


while the Rusanov equation becomes 


du , cdu 
dt dx 


- £_ cL S [ 5w - 4 - 15v 2 + + ••• (23) 


Kutler, bomax and Warming observed that the lowest order term 1. the 
modified equation representing dispersion or phase error was 
darivativ. term in the second-order case while in the third-order method, 
the lowest order diepersiv. term 1. a fifth derlvetive. For both method, 
the lowest order even, or dleslpetive. term is a fourth derivative. On. 
might expect that Improve.**, in -suing to the third-order technique 
should result fro. e decrees, in phase error. This conclusion was 
verified by Kutler. bomex end Weeing (DW tor .hocked flows. 
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Application to the case of non-shocked flows is discussed in a later 
section. 


The Kutler-Lomax-Warming Method Kutler, Lomax and Wanning (1972) 
have developed a noncentered version of the Rusanov scheme. There are 
two major differences between their technique and the original third- 
order Rusanov method. 

The K.-L-W method uses noncentered differences and uses the MacCormack 
method (evaluated at 2/3. AT ) for the first two steps while the Rusanov 
third level is used. This results in less computer storage and compu- 
tation time since only values at each mesh point are required. The 
other difference is that the fourth-order oo term has been differenced 


in a conservative manner so that the value of co can be altered during 
the calculations. If the co term is differenced as in the Rusanov method 
and altered during computation, incorrect wave speeds occur in the numer 
leal solution. 

The K-L-W method applied to the wave Eq. (G) takes the form 


,( 2 ) 

'J 


E 


n+1 


- f g [< l - £)F i +1 ■ < l - 2£)F " - «j-J 
i[ E j + E ^] - MS Hi! + (1 • 2 £)F j l) + <£ 

; - * J 


- l)F 


^ ft 


J+2 * 3E J + 1 + 3E ' 


a±m [ E " + i - 3E j + 3E "-i - E U] 


24 


- k M [- 2F "« + 7F ih - 7F "-i + 2 f j- 2 ] 

- 1 fi KS • F i-0 • 




( 24 ) 


Vl 


! 

i 

$ 

i 

4 

i 

t 

■i 


The form of this method applied to the two dimensional sample problem 
reads : 


,(D . A n 


i,j 

1< 2 > - - 
i,j 2 






- (vA) n . 

1>J_ 


,n+l 


i( A " j + |U«> 

2 l l, j i,j 3 lix [' i,j 

(w ) n 

-JLllIZla r n 

24 l A i+2,j 


(1 > - < u a) 5 1> 1 .1 -if [<«><» 

J i-l.JJ 3 Ay L l,.t 


,,n , _.n 4 n 3 

J 1+1 »J i,J 


A n 

( w x ) n 

+ [ A ?H • ’ 3A n . + 3A n . - A n .] 

24 l l+1 »J i.J i-l, J i-2 , j J 

(co ) n 

l.H/2 Ln . 3A n n n ] 

24 l.i > j + 2 3A i,j+l + 3A i,j A i,j.iJ 

* n - A n 1 
i,j-2j 


/ \1» 

(OJ } 

4^- i - - - 1 / 2 | A n . u n + 3A n 

24 L i,j+l 3A i,J 3A i,j-l 


_ 1 _ ( at r 

24 \ Ax |_ 


2( uA) 


i+2,j 


+ 7 (uA) 


n 


1+1. J 


- 7(uA)“ . + 2(uA) n .] 
i» J i-2, jJ 


+ % [- 2<v * ) ;, J « + 7<va >;, j+ i - ... + 2<» A > 


> n 

t.j-i 




3 /AT r, . . (2) 

• 8 (s L (uA) 


• <uA) i-Lj] 4 “ M!!* - 


i+i. j 


(25) 


where functionally 


H >n - w x k> l 

x i+l/2, J X L X 1+1/2, jJ 

, m %\ iv y ) 1 

y 1+1/2, J y L y 1 , J+1/2J 


( 26 ) 


( 1 ) 

i J- 
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and 


<,/ . 

x i+1/ 


= K (x >n + ( v n + ( v n , + Ux) i-i j I 

/2,j ‘L *1+2,3 i+l ’J i>J ,JJ 

n - ilix >n + <x v >n + C V" • + <Xy> t i-il 

i, J+l/2 U y 1,J« y 1,1+1 1,J lJ 


AT 

Ay 


( 27 ) 

This method ha. the .ame stability bound, and magnitude a. the 
Rusanov method and phase angle plot, shown in Figs. 3 and 4 for Rusanov'. 

method are applicable. 

Kutler, Lomax and Warming observed that this third-order method 
usually produced better results in problems involving large change, in 
the dependent variables than did second-order techniques. However, 
some sensitivity ol the results to variations In the local eigenvalue 
structure led them to attempt to "tune’' the method at each step by 

changing the value of the oi parameter. 

Application of the K-L-W method to the simple wave Eq. (10) result. 

in the modified differential Eg. (23) coamon to third-order Runge-Kutt. 
methods. To reduce the dissipation of the method, the to parameter should 
be selected to minimise the fourth derivative term. If the coefficient 
of this term is set equal to zero then 


/ 2 ,A 

(!) m U\f * V 


(28) 


which 1. lower stability bound a. predicted by linear theory. If 
minimum dispersion is desired, the coefficient o, the fifth-order ter. 

is set equal to zero giving 

(Ay 2 + 1)(4 - i/ 2 ). 

<0 • 3 


( 29 ) 
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Results presented by Kutier, Lomax and Wanning (1972) show that an im- 
provement in the solution is achieved if oo is altered according to 
either of these schemes. It should be noted that the value of v appearing 
in (28) and (29) is the local value of the Courant number. The and 
00y of (26) are thus determined through (28) and (29) with the appropriate 
values of v given in (27). 

The amplitude and phase error contour plots for the variable u) cases 
are presented in Figs. 5 through 8. The minimum dissipation amplitude 
plot in Fig. 5 shows that the "tuned" third-order method is less dissi- 
pative than the second-order technique. The effect of changing co to 
eliminate the first even term in the modified equation clearly reduces 
the damping. This should be compared with Fig. 7 which is also an ampli- 
tude contour plot but for minimum phase error as provided with to deter- 
mined through (26). It is apparent that the damping is greater in this 
case but only below Courant numbers of about 0.6. This should have 
been anticipated since the allowable range of w is restricted at higher 
Courant numbers and the method becomes less sensitive to the value of to 
used. 

The phase angle contours show that the minimum dispersion case in 
Fig. 8 provides a smaller phase error over the entire range of stable 
values of v than the minimum dissipation case shown in Fig. 6. On the 
basis ox the amplitude and phase angle plots one would expect results 
using the "tuned" K-L-W method to be superior to the constant to third- 
order methods. Application to a problem justifying this conclusion is 
presented in the following section. 


3. Numerical Calculations 

A simple velocity field is prescribed over a coordinate mesh if 
it is assumed that it is produced by solid body rotation about some axis. 

In this case the divergence of the velocity field is zero and the con- 
servative form given in (3) is valid. The initial distribution of A is 
chosen to be a cone of unit altitude at the center decreasing to zero 
at the base radius of four cell widths. The following data were held 
constant during the calculations unless otherwise noted: 

n * number of time steps * 40 

(7* rigid body angular velocity = - .001 radians/sec 
At ■ time step » 30 seconds 
Ax-Ay«A*1.0* mesh size 

Center of rotation (12, 12) 

Center of initial A distribution (18, 12). 

A schematic representation of this problem is shown in Fig. 9. 

Figure 10 presents the results of using MacCormack's method with 
a forward predictor, backward corrector while the results obtained with 
a backward predictor, forward corrector are shown in Fig. 11. The results 
in both cases are similar in terms of correctly advecting the center of 
the disturbance and in reproducing contours of constant A. The results 
are very similar to those produced using the Lax-Wendroff method as 
reported by Fattahi (1971). 

The Rusanov third-order method using to * 3 was applied to this prob- 
lem and the results are shown in Fig. 12. These results are not very 
encouraging. However, a more detailed analysis of the situation produces 
soma important information. The Rusar.jv method has been shown to produce 
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results which are unfavorable when the Courarit number is not close to 
one and the upper limit value of 3 is used for co- In this problem the 
mesh ratio is AT/A * 30 which gives an effective Courant number of 
approximately 0.35 based on the maximum velocity. This means that the 
local Courant number at each point in the mesh is less than or equal to 
0.35. Referring to Figs. 3 and 4, the i.'.ost dissipative and dispersive 
regions of the Rusanov method lie in this range. This is the reason for 
the poor results shown in Fig. 12. It is interesting to note that the 
MacCormack method provides a better solution for this case than the 
third-order Rusanov technique. Figs. 1 and 2 show that less dissipation 
and dispersion are introduced in the solution using the MacCormack method 
for the range of effective Courant numbers encountered in this numerical 
example. For this case the minimum value of co for stability required 
by the Rusanov method is approximately 0.5 while the maximum value is 3. 
Calculations using values of c0 near the lower bound have resulted in 
good solutions for problems involving shocked flows like those of Kutler, 
Lomax and Wanning (1972) «uu also for the sample problem presented here 
by Burstein and Mirin ( 1970 ) and Fattahi. Values of cd near or below 
the linear stability bound result in excellent solutions for shocked 
flows. Fattahi ( 1971 ) reports a similar conclusion and, in solving 
the exact problem as presented here, notes that an cd value of about 0.6 
produces the best results. This is shown in Fig. 13 . Solution* obtained 
with uo below the linear stability bound provide excellent results for 
short integration times. For long integration times the solution 
eventually becomes unstable if a value of cc outside the allowed range 


is used. 
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The problem of selecting the proper value of tots alleviate! 
the Kutler-Lomax -Warming technique with either minimum dispersion or 
dissipation 1. used. Figures U and 15 present the results of these 
calculation, for a time step of AT - 30 seconds. The minimum dissipation 
case provide, better results at this value of AT. Referring to the am- 
plitude contours for these two methods reveals that the damping is 
severe for the minimum dispersion case at low Courant numbers. This is 
apparent if the maximum adverted value is compered for these two cases. 

This type of problem is evidently more sensitive to dissipative 
than error, in phase angle. Kutler, Lomax and Warming (1972) concluded 
that the minimum dispersion provided better results for their work 
involving shocked flows. The difference in the results using the two ways 

of calculating a diminishes as the linear stability bound is approache 

• ► n nf AT = 75 sec are shown in Figs. 16 and 17. As 

Solutions for a time step of AT - /h sec are 

noted from the damping and phase error curves, the results obtained should 
be similar. However, the minimum dissipation solution still appears to 


be more satisfactory. 


4. Conclusions 


A comparison of relative errors and computation time, are presented 
in Table 1. Based on accuracy in terms of maximum adverted value, phase 
error and error in the r.dlu. of the center location, the third-order 
method, sr. usually better unless relatively low value, of Courant 
number are coupled with large .slues of » The third-order method, 
require the .... magnitude of computer time a. HacCormack's method. 

They ere still an order of magnitude faster than the fourth-order space. 


«,*-<** TM ffcMl tig 
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second-order time Roberts and Weiss method. Although this technique was 
not investigated here, Molenknmp has reported the re La tive time 

required for its use. 

Of the methods tested, the Kut Ler-Lomax-Warming method tuned for 
minimum dissipation provides the best comparison with the exact solu- 
tion. This method should be studied in more detail. In particular 
results using the various combinations of the c's defined in (14) 
should be compared. The results reported here used a forward predictor 
and backward corrector to advance to 2/3AT and central differencing in 
the third step. In addition, application should be made to more com- 
plex problems to assess the usefulness of the technique of altering 
the stability parameter cj during the computation. While the eigen- 
value structure required in choosing the proper t C in a complex problem 
is difficult, reasonable approximations can be made. The increase in 
accuracy of the resulting solution far outweighs the additional effort. 
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Table 1. Comparison of Calculated Results 
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Fig. 1. Constant |5| contour for MacCormack ( s second order method. 

Fig. 2. Contour of - 0/0a for MacCormack's second order method. 

Fig. 3. Constant \$\ contours for Rusanov Method with 0) * 3.0. 

Fig. 4. Phase angle for Rusanov r.jfchod with <*> * 3.0. 

Fig. 5. Constant |§| contours for K-L-W Method for minimum dissipation. 

Fig. 6. Relative phase error for the K-L-W Method tuned for minimum 
dissipation. 

Fig. 7. Constant j?| contours for the K-L-W Method with minimum 
dispersion. 

Fig. 8. Contours of - 0/'0a for K-L-W Method with minimum dispersion. 

Fig. 9. Schematic representation of the advection problem showing 
the initial A distribution. 

Fig. 10. Solution using the MacCormack (forward predictor, backward 
corrector) Method • 

Fig. 11. Solution using the MacCormack (forward predictor, backward 
correc tor) Method • 

Fig. 12. Solution using Rusanov's Method with w = 3.0. 

Fig. 13. Solution using the Rusanov-Burstein-Mirin Method, u) - 0.6, 

AT « 30. 

Fig. 14. Solution using the K-L-W Method tuned for minimum dissipation 
AT * 30 sec. 

Fig. 15. Solution using the K-L-W Method tuned for minimum dispersion 
AT ■ 30 sec. 

Fig. 16. Solution using the K-L-W Method tuned for minimum dissipation 

AT - 75 sec. 

Fig. 17. Solution using the K-L-W Method tuned for minimum dispersion 
AT - 75 sec. 
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